Multiresolution analysis for efRcient, high precision all-electron density- functional 

calculations 
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Multiresolution analysis of electronic structure affords the opportunity to capture the full physics 
of atomic cores in a systematically improvable manner. Applying new techniques, we demon- 
strate for the first time that multiresolution analysis of all-electron calculations within density- 
functional theory can be carried out to high precision with a computational effort comparable to 
that of the corresponding plane-wave pseudopotential calculation, which neither captures the full core 
physics nor is systematically improvable. With this approach, we present calculations of paramag- 
netic core-level shifts where local density-functional theory is the sole uncontrolled approximation. 
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I. INTRODUCTION 

Over the last several decades, the ab initio density- 
functional approach has proven an accurate, reliable and 
effective tool for the study of condensed matter. It has 
found application in such diverse areas as the study of 
surfaces, point defects, melting, diffusion, plastic defor- 
mation, disorder, catalysis, phase transitions and chem- 
ical reactions H. In principle, the only approximation 
required in the practice of density-functional theory is 
some model for exchange and correlation effects. The fact 
that the forms for exchange and correlation are universal 
and independent of a priori knowledge of the physics or 
chemistry of the systems under study gives reason for far 
greater confidence in the predictions of density- functional 
calculations than those of their semi-empirical counter- 
parts. 

Present practice of electronic structure calculation, 
however, falls short of this ideal. All approaches used 
in current production employ prior knowledge of chem- 
istry as a criterion for freezing out degrees of freedom in 
order to deal with the special demands of representing 
physics near the Coulomb singularity of the atomic nu- 
cleus. While such schemes succeed in making ab initio 
studies of complex phenomena accessible, they make it 
difficult to improve the calculations systematically and 
face the danger of introducing artificial biases. For in- 
stance, Gaussian basis sets are built directly from chemi- 
cal intuition of how atomic orbitals are expected to polar- 
ize g. The hncar muffin tin orbital (LMTO) method g, 
the linearized augmented plane wave (LAPW) method 
@, and the full potential LAPW (FLAPW) method |] 
all use one type of basis set inside of a set of spheres 
organized around the nuclei and another type of basis 
set outside of the spheres, thus treating physics differ- 
ently in different regions. Finally, the pseudopotential, or 
"frozen core", approximation, ignores potentially impor- 
tant effects such as core polarization and interferes with 
the valence wave functions. As a result, it is not uncom- 
mon for different groups using different techniques, or 



even the same technique, to disagree on the predictions 
of density-functional theory. 

If one could perform ab initio calculations with compa- 
rable computational effort, but without the biases of the 
above approaches and with the ability to systematically 
improve the basis, one finally could access directly the 
predictions of density- functional theory without ambigu- 
ity. Also, one could better explore the development of im- 
proved energy functionals, which demands not only bet- 
ter functional but also the use of basis sets with highly 
controlled precision sufficient to resolve reliably the asso- 
ciated millihartree- level improvements. 

Indeed, wavelet bases do not incorporate any prior 
knowledge of electron physics and provide a natural 
framework for systematic expansion of the high spatial 
frequencies which electrons exhibit near nuclei. Pio- 
neering all-electron density-functional calculations based 
on multiresolution analysis have been reported, first on 
atoms Q and then molecules [0J8|. Other applications 
of multiresolution analysis to related problems include 
single-electron problems |9| , pseudopotential calculations 
pO|-p^ and the solution of Poisson's equation Q . The 
initial applications to electronic structure, however, were 
based on relatively primitive techniques and were never 
demonstrated to produce results at the limits of the ac- 
curacy of density-functional theory nor to provide results 
with an efficiency comparable to standard techniques 
such as the pseudopotential plane-wave approach. We 
believe that the reason for this shortcoming is that stan- 
dard approaches in the wavelet literature are designed for 
applications, such as digital signal processing, which have 
far different demands than do continuum problems in the 
physical sciences. In response, we have taken a different 
tact. Starting with the fundamental principle underly- 
ing wavelet theory, multiresolution analysis |l^-fi6|, we 
have developed a unique set of methods and tools de- 
signed specifically for the treatment of continuum prob- 
lems P,kIJi7|- Below, we report the results of the use of 
these new techniques. 



II. CHOICE OF BASES 

Figure |l| illustrates the general structure of the mul- 
tiresolution analyses employed in this work. These begin 
with a coarse representation consisting of a uniform or- 
thorhombic grid of rather extended basis functions (large 
circles in the figure) placed in a unit cell with periodic 
boundary conditions (dashed square). Additional basis 
functions of intermediate extent (medium circles) on a 
grid of one-half the original spacing then carry details 
which the coarse representation does not. Yet smaller 
functions (small circles) on a grid of one-half the spacing 
of the intermediate grid carry information for one further 
scale down. Repetition of this construction provides for 
any desired level of resolution. 

With appropriate choice of basis functions, multireso- 
lution analysis ensures the crucial result which allows the 
following calculations to achieve systematic convergence, 
that such a multiscale basis spans exactly the same func- 
tion space as would a uniform basis of functions on the 
finest scale. The basic idea is to insist that each basis 
function be expressible exactly as a linear combination of 
translates of itself scaled down by a factor of two, 



necessary only in the region closest to the nucleus (small 
solid square), and the intermediate scale functions are 
necessary only in a region extending somewhat further 
from the nucleus (medium solid square) . Restricting the 
multiresolution analysis by maintaining only those func- 
tions (filled circles) within the refinement regions for each 
scale affords significant savings while affecting the fi nal 
results negligibly — shortly below in Section IIIC| , we 
demonstrate the effects of such restriction to be control- 
lable down to at least the TTiicrohartree per atom level. 

To describe fully the bases used in our calculations, 
we must specify both the original multiresolution analy- 
sis and its restriction. All multiresolution analyses em- 
ployed in this work begin with a coarsest level of one 
bohr spacing and employ basis functions of the "third- 
order semicardinal type" . This means that the functions 
satisfy the two-scale relation (|l|) with the two additional 
constraints of "cardinality," 



Hn) = ^s,o 



(2) 



b{x) = y^ Cfib{2x — n) , 



(1) 



(zero value on all integer points but n = 0), and "exact 
reconstruction" of multinomials up to third order, 
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n"nf^n2b{x ~ fi); 0<a,/3,7,<3. (3) 



where b(x) is the basis function and n ranges over all 
triplets of integers. The above equation is known as the 
two-scale relation. Applied recursively, this relation im- 
plies that all coarser functions in the basis are expressible 
exactly as linear combinations of functions on the finest 
scale and, therefore, that the multiscale basis is equiva- 
lent to the basis on the finest scale. For a full review, see 
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FIG. 1. Schematic representation of a two-dimensional re- 
stricted multiresolution analysis of three levels of resolution. 
(See text for description.) 

The basis consisting of all functions indicated in the 
figure contains precisely as many functions as the equiva- 
lent uniform basis on the finest scale. Thus, its direct use 
would represent no savings. However, the coefficients of 
the basis functions on the finest levels of the multiresolu- 
tion basis represent high frequency details and hence are 
significant only in the vicinity of the nucleus (diamond in 
the figure). Consequently, the finest scale functions are 



These two additional constraints, respectively, simplify 
the computational algorithms considerably and allow 
third-order intcrpolative representation of functions such 
as the electron density. See [0| for a full discussion. 

To describe the restrictions employed in this work, we 
note that they consist of orthorhombic refinement re- 
gions. Thus, specification of the restriction reduces to 
identification of the dimension {Nx,Ny,Nz) and location 
of these regions at each level of resolution. Because the 
refinement regions of a finer level (higher level number) 
are always contained within the refinement region of the 
next higher level, it is most convenient to identify the lo- 
cation of a refinement region in terms of the displacement 
of its origin from the origin of its parent region in units of 
the parent's spacing (A^;, A^, A^). Using this approach, 
for example. Table J specifies the two-dimensional restric- 
tion which Figure [I illustrates. 

In some cases, such as the restriction for the calcu- 



lation of the oxygen molecule in Table IV, a single re- 
finement region contains multiple children, regions of the 
next level of refinement. In such cases, we first list the 
first child with all of its progeny directly underneath and 
then proceed to the remaining children of the original 
parent. Ordered and listed with the data as appear in 
these tables, such a specification is sufficiently general, 
complete and compact that we use it directly as an input 
data structure when specifying multiresolution analyses 
to our software. 



III. EFFICACY OF MULTIRESOLUTION 
ANALYSIS OF ELECTRONIC STRUCTURE 

A. Kohn-Sham equations in a multiresolution basis 

For a full review of the form and solution of the Kohn- 
Sham equations in a multiresolution basis, see ||7|]. Here, 
we sketch briefly the key equations which we use to pro- 
duce the results in this work. The central quantity in 
density-functional theory is the energy expressed in terms 
of the expansion coefficients Cai of the Kohn-Sham or- 
bit als. 



■0i(x) 
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(4) 



where i runs over the Kohn-Sham orbitals and a ranges 
over the functions in the basis employed in the calcula- 
tion. In [0JlJ] , we give the appropriate expression for this 
energy function in an arbitrary basis. This expression in- 
volves standard matrix operations, such as the trace Tr , 
and the basis-dependent matrices 



Oci3= / ba{x)*bi3{x)£'x 



(5) 



Lai3= / ha{x)*V'^hp{x)d^x 

where a and (3 range over all functions in the basis and 
p ranges over a set of sample points in real space, which 
in the present case are the centers of the wavelet basis 
functions (filled circles in Figure nh. The first two matri- 
ces give the overlap integrals and matrix elements of the 
Laplacian, respectively. The final matrix consists of the 
values of each basis function at each sample grid point p 
so that the product of X with a vector of expansion coef- 
ficients transforms it to a vector of values in real space. 
In terms of the above matrices, the density-functional 
expression for the total energy in terms of the expansion 
coefficients gathered into the matrix C is 



Elda{C) = Tr (fCH-^L)c) + (I-^n)^Ov 



(6) 



Here, / is the occupancy of each orbital (typically 2), u is 
the vector of expansion coefficients of the ionic potential. 



V,c 



y^^Vaba{x), 



£xcin) is the usual exchange-correlation energy of a ho- 
mogeneous electron of uniform density n, and n and d, re- 
spectively, are the electron density evaluated at the sam- 
ple points and the expansion coefficients of the Hartree 
potential. 



n = diag((JC)/(JC)^) 



(7) 
(8) 



where "diag" is the operation of forming a vector from 
the diagonal elements of a matrix. The final quantity 
needed to solve the Kohn-Sham equations is the deriva- 
tive of Elda with respect to the orbital coefficients C, 



dE 



dC 
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(9) 



{j-''Od)}lC, 



where Diag a represents formation of a diagonal matrix 
with diagonal elements given by the components of the 
vector a. 

To aid the reader's interpretation of (|6|-g), we note 
that these expressions are fully general and applicable 
not only to wavelets but also to plane waves. These 
expressions, therefore, represent minor generalizations 
of the same sequence of operations found typically in 
plane-wave density-functional codes. For instance, in the 
computation of the real-space electron density (0), the 
quantity {IC) represents the fast Fourier transformation 
of the orbitals from coefficient to real space, the outer- 
product (XC) f {XCy is then the real-space density ma- 
trix, whose diagonal elements ultimately give the real- 
space electron density n.R Similarly, to find the Hartree 
potential (||), one first multiplies X^^n (inverse Fourier 
transforming the real-space density to coefficient space in 
the plane wave case), multiplies by O (for plane waves, a 
simple volume normalization factor) and then by —An 
times the inverse of the Laplacian matrix (47r/G^ for 
plane waves). The terms in Eqs. (0,H) each have simi- 
lar interpretations. 

To produce the results reported below, we minimized 
the expression (^) over all possible coefficients C, using 
the diagonally preconditioned conjugate-gradient algo- 
rithm augmented with the analytically continued func- 
tional approach pii| to handle the orthonormality con- 
straints on the orbitals. Combined with (j^-H), this de- 
scribes our calculations fully and explicitly. 

None of the calculations below employ gradient cor- 
rected density functionals. During the preparation of the 
manuscript, however, we were asked to comment on how 
one would handle the numerical issues which such func- 
tionals raise. Our general approach to numerical issues 
is to note that so long as the evaluation of a given term 



^Eq. (|7|) is a formal expression. Anticipating the diag oper- 
ator, practical software computes only the diagonal elements 
and not the full density-matrix. Similarly, in Eq. IM), rather 
than forming large matrices such as Diag a explicitly, our soft- 
ware computes products such as (Diag a)b by multiplying each 
element of b by the corresponding element of a. 



in Eloa is exact for any finite expansion of the Kohn- 
Sham orbitals (jj), then all numerical issues reduce to the 
quality of the expansion for the wave functions, which we 
show below to be well-controlled in multiresolution bases. 
Therefore, to avoid numerical instabilities associated 
with finite differencing to evaluate quantities such as 
Vn{x), we would evaluate the exact, analytic gradient 
of the charge density associated with the orbital expan- 
sion (H). To do this, we define an additional matrix 

Gpa = ^ba{p), 

where a and p range as in the definitions (g|). The values 
of V'0i(a') at the sample points p are then 



vMp) = [QC] 



pi '■ 



so that the gradient of the charge density at these points 
is exactly 



i 
i 

= [(jc)/(gc)t-f(gc)/(JC)t]^^ 

Vn = diag ((JC)t/(eC) + (gC)V(IC)) • 



When the exchange-correlation function is evaluated us- 
ing this quantity and inverse-transformed with I~^, the 
resulting coefficients, as per the discussion in Section |IV|, 
will be exactly the same as would be obtained were 
the numerical evaluation carried out in a full multires- 
olution analysis at arbitrary resolution without restric- 
tion, thereby mitigating any numerical issues associated 
with the evalatuation of gradient corrected functionals 
on wave functions expanded as in (0). 



B. Systematic convergence 

Present feasible approaches to all-electron calculations 
complicate access to the full predictive power of density- 
functional theory and the development of new function- 
als because these methods approach the atomic cores and 
the valence regions in fundamentally different ways and 
therefore are difficult to bring to convergence in a sys- 
tematic manner. While it has been clear for some time 
that the ability of multiresolution analysis to focus reso- 
lution in small regions of space can be exploited to per- 
form all-electron calculations |^-gl , the highly systematic 
nature of multiresolution analysis and the consequent ad- 
vantages have yet to be explored in the context of elec- 
tronic structure. We now present the first demonstra- 
tion of a multiresolution analysis reproducing all-electron 



density-functional results with highly systematic conver- 
gence down to millihartree accuracy. 

For this demonstration, we compare results for atoms 
as calculated within the local-density approximation as 
parameterized in pO| using both standard techniques and 
our multiresolution approach. Atoms provide a conve- 
nient test bed because spherical symmetry produces an 
effective one-dimensional problem which standard tech- 
niques then readily solve to a precision controllable to 
better than 1 nanohartree. To provide a fair measure 
of how we expect the multiresolution method to perform 
in practice, all calculations with multiresolution analyses 
in this work are fully three-dimensional and without any 
simplifications. 

To underscore that our approach is feasible beyond first 
and second row elements, we have chosen calcium for 
comparison. To represent the charge distribution of the 
calcium nucleus (and all nuclei in this work), we use the 
linear combination of three Gaussian distributions with 
the widths and charge contents appearing in Table 0. 
This charge distribution reproduces the atomic electron 
eigenvalues to within 0.5 mH [ pl| , an accuracy which 
is systematically improvable through a renormalization 
approach [p2|. All calculations and comparisons below 
are made for this form of nuclear charge distribution. 
As described above, the multiresolution analysis begins 
with a spacing of 1 bohr on level zero. To explore the 
convergence toward the final solution, we perform self- 
consistent calculations with from seven to nine levels of 
refinement, for a resolution of down to 1/2^ ~ 0.002 bohr. 
Table III summarizes the restrictions employed for the 



calculation. 



Calcium: total energy [hartree] 
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FIG. 2. Error in total all-electron energy within the local 
density approximation, A_E, of a multiresolution calculation 
of an isolated calcium atom at 7, 8 and 9 levels of resolution 
(squares) below a top level of 1 bohr spacing. The horizontal 
axis gives the cube of the resolution on the finest scale (Ax)"^. 

Figure Q shows, as a function of the cube of the spacing 
on the finest level, the difference between the total en- 
ergy calculated with the multiresolution analysis and the 



highly precise solution of the equivalent one-dimensional 
problem. When we employ nine levels of refinement, the 
absolute error is 6 mH, one part in 10^ of the total energy. 
The simple, linear approach of the error toward zero 
in this plot is a direct consequence of the fundamental 
nature of the basis. The multiresolution analysis which 
we employ satisfies the two-scale relation and, thus, is 
equivalent to a uniform basis on the finest scale of reso- 
lution by construction. The linear behavior evident in the 
figure, therefore, can be understood a priori as a direct 
consequence of the third-order nature of the semicardinal 
basis. The only other method which provides this level of 
analytic understanding of the errors as a function of the 
size of the calculation is the plane- wave approach. How- 
ever, large scale plane- wave calculations are only feasible 
in practice with pseudopotentials and therefore do not 
systematically converge to the correct all-electron result. 
Multiresolution analysis, therefore, is unique among fea- 
sible all-electron approaches in allowing for a high degree 
of a priori understanding of and systematic control over 



C. Uniformity of description of space 



function from the coarsest scale and then move the atom 
along a coordinate axis in increments so that the nucleus 
falls directly on the centers of basis functions, eventually 
sampling functions from each of the eight levels in the 
calculation. As the atom moves, we add or remove func- 
tions from the basis to maintain the atom at the center 
of each refinement grid while conserving the number of 
basis functions on each level. 

Figure shows the results of the above calculations. 
Despite the switching on and off of the basis functions 
and the radically changing nature of the basis function 
on which the nucleus falls (ranging by over two orders of 
magnitude in length scale), the root-mean-square fluctu- 
ation in the energy is less than 3 microhartree, which we 
could reduce yet further by expanding the extent of the 
refinement regions. In terms of Pulay- force corrections 
PSJ , the maximum fluctuation in energy as we move the 
atom by 1/64 bohr (the smallest step for which the basis 
actually changes) is 2.0/iH, corresponding to a maximum 
average Pulay force of ~ 128 ^H/B ?« 0.007 eV/A. This 
remarkable result makes the multiresolution approach an 
extremely attractive tool in calculations where the atoms 
move, such as structural relaxation or molecular dynam- 
ics. 



Present methods demonstrated to be viable for all- 
electron calculations all treat space near the atomic nu- 
clei differently from the interatomic regions. This inher- 
ently biases the description toward the location of the 
atoms and complicates the evaluation of forces on the 
atoms [E3|. Surprisingly, even under restriction, this ef- 
fect is nearly absent in multiresolution analysis because 
of the unique way in which multiresolution analysis rep- 
resents information. 

In direct contrast to the expansion coefficients in finite- 
element bases or finite-difference calculations, which are 
in proportion to the value of the represented function in 
the corresponding region of space, the expansion coef- 
ficients of the finer-scale functions in a multiresolution 
analysis represent only the higher-frequency details in 
the corresponding region. The expansion coefficients of 
the finer-scale functions which are restricted from a mul- 
tiresolution basis (empty circles in Figure |l|) , unlike the 
coefficients which would appear in those regions in tradi- 
tional approaches, therefore can be made to vanish con- 
troUably by expanding the corresponding refinement re- 
gions (squares in Figure ^ outward from the nuclei. As 
the coefficients dropped from the expansion vanish, the 
restricted multiresolution analysis becomes indistinguish- 
able from the full analysis, and the underlying description 
of space therefore becomes uniform. 

To demonstrate this, we consider the calculation of the 
total energy of a single atom (aluminum) at different lo- 
cations in space. For these calculations, we employ a 
cubic cell of dimension 48 bohr with seven levels of re- 
finement, each consisting of a cubic array of 48'^ basis 
functions. We begin with the atom centered on a basis 
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FIG. 3. Fluctuation {AE) in total all-electron energy of an 
isolated aluminum atom as a function of distance (Ax) of the 
nucleus from the center of a basis function on the coarsest 
scale. Note that vertical scale is in microhartree. 



D. Description of chemical bonding 

To investigate the efficacy of the multiresolution ap- 
proach in the description of chemical bonding, we con- 
sider the oxygen molecule. Given the paramagnetic na- 
ture of this molecule, we carry out the calculation within 
the local spin-density approximation (LSDA) as param- 
eterized in [g^. To ensure sufficient isolation of the 
molecule from its periodic images, we employ a super- 
cell of dimensions 14.8Axl2.7Axl2.7A. We then orient 
the molecule along the long axis of the cell, place the 



nuclei at their experimental separation, and employ the 
multiresolution analysis in Table IV . 

For comparison, we also performed, within the same 
geometry and parameterization of the local spin-density 
functional, a plane-wave calculation within a periodic 
supercell of identical dimensions using a pseudopoten- 
tial optimized for convergence according to the proce- 
dure of Rappe et al. ]25[| . This pseudopotential required 
a plane wave cut-off of 35 hartree to converge the to- 
tal pseudo-energy to within 0.10 hartree. Given this 
level of convergence in the total energy and the fact that 
the fractional error in the variational total energy scales 
like the square of the fractional errors in non-variational 
quantities such as the eigenvalues, we estimate the in- 
dividual eigen-energies to be converged to within about 
0.03 hartree at this plane-wave cut off. 

Figure ^ compares the Kohn-Sham eigenvalues from 
the pseudopotential calculation with those from multires- 
olution analysis. Within the expected convergence of the 
pseudopotential calculation, the agreement is complete. 
These results illustrate that the multiresolution approach 
gives a description of bonding at least as reliable as the 
standard pseudopotential approach. Because the issue of 
transferability remains an open question in pseudopoten- 
tial theory, this result also lends greater credence to the 
pseudopotential which we have constructed for these cal- 
culations. We envisage the use of multiresolution analysis 
in this way to aide in the construction of pseudopoten- 
tials with greater transferability. 



Energy (Hartree) 
i MRA 



-0.4 



-0.6 



-0.8 



-1.0 



-1.2 



PW 



AE=0.002H 



AE=0.020H 
AE=0.009H 



AE=0.032H 



AE=0.010H 



FIG. 4. Kohn-Sham eigenvalues for valence states of the 
oxygen molecule: results of all-electron multiresolution anal- 
ysis calculation (MRA) and plane-wave pseudopotential cal- 
culation (PW). 



IV. COMPUTATIONAL EFFICACY 

The computational demands of wavelet-based all- 
electron calculations in previous works were so extreme 
that the literature has yet to explore the practicality of 
the computations. The calculations presented in this 
work are the first to apply a series of new techniques 
[f7|, p7|j26| to achieve high performance in all-electron prob- 
lems, placing us in a position to account in detail for 
the computational effort that our calculations demand 
and thereby provide a standard against which the per- 
formance of future wavelet calculations can be compared. 
We find that our new techniques make the demands of 
all-electron calculations sufficiently comparable to those 
of frozen-core pseudopotential calculations of the same 
precision to make wavelet calculations a practical alter- 
native. 

Rather than using techniques developed for other 
classes of application, the techniques which we use here 
are tailored to the specific demands of continuum prob- 
lems in the physical sciences. Our basis, for instance, 
technically is not a wavelet basis because its dual con- 
sists of sets of Dirac-delta functions. These delta func- 
tions, however, are precisely what are needed to recover 
exact representations from samples of physical fields at 
extremely limited numbers of points in real space p7[ . 
In the present context, this means that the operations 
I^^ appearing in (^,||J9|) always yield precisely the same 
results as would be obtained were the calculation carried 
out with a full, unrestricted multiresolution analysis of 
arbitrary resolution. (See Q,|l3l-) ^s evident in (g,||,||), 
such operations are key in the effective treatment of non- 
linear couplings such as exchange and correlation effects 
and evaluation of the Hartree potential. Without this 
profound result, accurate evaluation of such eflfects re- 
quires evaluation of quantities at a much larger set of 
points in real space, thereby slowing the calculation con- 
siderably. To go along with these new methods for the 
transforms {I, X~^ and related operations), we have also 
developed a class of 0{N) methods for multiplication by 
the operators O and L with similar exactness properties 
[l7|,[l7| . Finally, to further improve the performance of the 
calculations, we have developed special cache-optimized 
algorithms for all of these methods ||2^ . 

As a case study, we consider in detail the computa- 
tional demands of the calculations of the electronic struc- 



ture of the oxygen molecule presented in Section HID 
above. To make the comparison with a pseudopoten- 
tial calculation meaningful, we choose a restriction of the 
multiresolution analysis (Table E^) and a plane- wave cut 
off (35 H) which achieve similar convergence in the total 
energy: 0.08 H and 0.10 H, respectively. 
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FIG. 5. Convergence of total energy of the oxygen molecule 
as a function of actual wall-clock time on a 400 MHz Pen- 
tium II PC for basis sets reaching a similar level of con- 
vergence: plane-wave pseudopotential calculation (dashed 
curve), all-electron multiresolution analysis calculation (solid 
curve) . 

Figure || presents the convergence of the total energy 
as a function of actual wall-clock time on a 400 MHz 
Pentium II PC both for our wavelet calculation and for a 
plane-wave pseudopotential calculation performed with 
the highly optimized DFT-|-f code ||l|]. Despite the 
fact that the pseudopotential calculation freezes out the 
physics of the core, the absolute convergence for the pseu- 
dopotential calculation is slower than that of the all- 
electron wavelet calculation. These results are the first 
to establish multiresolution analysis as a computationally 
viable, unbiased and systematic approach to the calcula- 
tion of all-electron systems. 

Several advances have taken place to make possible 
this favorable comparison. To place these in context, we 
consider the time required for the single most time con- 
suming basis dependent operation in the two calculations, 
the Laplacian for the wavelet case and the Fourier trans- 
form for the plane-wave case. In general, the total time 
required to perform such operations is 



T = Na 



NbNf/i 



N 



F/T 



(10) 



where Na is the number of times the corresponding op- 
eration is performed while reaching the desired level of 
convergence, Nb is the number of functions in the basis, 
Np/B is the number of floating point operations needed 
to process each basis function during the operation, and 
Np/rp is the number of floating point operations which 
the algorithm for the operation achieves per unit time. 

Table ^ presents these quantities for calculations 
which achieved the convergence of 10~® H in Figure ra. 
The table confirms the finding established in earlier works 
that the total number of basis functions Nb required for 
both calculations is similar ||,0l. 

The new algorithms which we employ Mu^ exploit 
the local nature of the multiresolution basis to produce 



a data fiow pattern which now is sufficiently simple to 
reduce Np/B from that of our previous calculations jg] 
by several orders of magnitude. To multiply by X, for 
instance, we apply the two-scale relation (|l|) recursively 
until we reach an expansion in terms of functions on the 
finest scale, which, by virtue of (||), gives the values of 
the function at each point p. Eq. (nl) implies that each 
step in this recursion is simply a three-dimensional con- 
volution by the sequence Cfj. Reference presents meth- 
ods for wavelet bases which similarly decompose multi- 
plication by all remaining basis dependent matrices (O, 
L, I~^ , l\ ^~'^) into three-dimensional convolutions for 
wavelet bases. For our choice of basis, these convolutions 
are short and separable, leading to an operation count 
with a significant advantage over that of the fast Fourier 
transform: Np/B ~ 250 versus A^p/^ « 751og2(15A^) (in- 
cluding appropriate factors for mapping the plane-wave 
sphere into the Fourier transform box). Table VI shows 



that for the present calculations, the wavelet Np/B is 
superior by a factor greater than six. 

The simpler communication pattern of separable three- 
dimensional convolutions also allows us to develop spe- 
cial algorithms with significantly improved cache perfor- 
mance and thus processing rates Np/p nearly an order 
of magnitude improved over that achieved in previously 
reported wavelet-based electronic structure calculations. 
Table W\ shows that, in fact, we now can achieve floating 
point operation rates Np/p somewhat superior to those 
of the highly tuned FFTW Fourier transform package 
PTj , which we employ in the plane wave calculations. Fi- 
nally, we note as a possible area for future research that 
the only area where the plane wave calculation holds a 
significant advantage is in the number of applications re- 
quired of the rate limiting operator Na- 

The result of all of the above factors is that the plane 
wave and wavelet calculations require comparable (within 
a factor of three) amounts of time T in their most time 
consuming basis-dependent operators. The remaining 
components of the calculations do consume significant 
amounts of time and so the total run times (Figure @) 
are significantly longer. In the wavelet case, our data 
for the Laplacian operator accounts for about one-half 
of the total run time. The additional operations in the 
plane wave calculations consume a much more significant 
fraction of the total time. As a result, in this particular 
case it turns out that the final total run time for the 
pseudopotential calculation, which does not include the 
physics of the core explicitly, is in fact noticeably longer 
than that for the corresponding wavelet calculation. 



V. CORE-LEVEL PHYSICS 

One great advantage of multiresolution all-electron cal- 
culations over their pseudopotential counterparts is the 
direct access which they give to the full physics of the 
atomic cores in an unbiased representation. A promising 



area of application for such calculations is their use as 
an aid to the interpretation of experimental techniques 
which measure environmentally dependent shifts in the 
core states, such as electron energy loss spectroscopy 
(EELS) [p8[ and electron spectroscopy for chemical anal- 
ysis (ESCA) Ig^. As a simple demonstration of the 
physics accessible through this approach, we consider the 
environmental dependence of the oxygen Is state in the 
water and the diatomic oxygen molecules. 
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FIG. 6. Difference, along the symmetry axis of the 
molecule, in total charge density of the oxygen Is states (both 
spin channels) in going from atomic oxygen to the oxygen 
molecule. Values expressed as fractions of the peak density of 
the Is state in atomic oxygen. Vertical lines indicate locations 
of the oxygen nuclei. 

For diatomic oxygen, which is paramagnetic, we use 
the results of the local s pin-density calculations described 
above in Section HID. For the water molecule and the 



reference atomic states, we employ the local density ap- 
proximation as parameterized in EQ]. For all calcula- 
tions, we place the nuclei at their experimentally known 
locations and employ the multiresolution analyses which 
Tables [V and M summarize. The general, simple struc- 
ture of multiresolution analysis allows this brief descrip- 
tion to specify the calculations completely. 

Our results for the deformation of the Is state of ox; 
gen in the above molecules appear in Figures and 
Figure g shows the change in the total electron density 
in the oxygen Is states (sum both over spin channels 
and over bonding and anti-bonding states) in going from 
atomic oxygen to the oxygen molecule. Despite being 
quite small (at the 0.03 % level), the core polarization 
effect is quite clearly defined in our multiresolution anal- 
ysis, which we emphasize builds no a priori knowledge 
of this polarization into the calculation. Interestingly, we 
find that, despite the buildup of repulsive valence charge 
in the bond between the oxygen atoms, the decreased 
screening of the nuclei is sufficient to produce local elec- 
tric fields which polarize the core electrons toward the 
center of the molecule. 
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FIG. 7. Difference in charge density of the oxygen Is state 
in going from atomic oxygen to the water molecule along a 
line directly toward one of the protons. Values expressed as 
fractions of the peak density of the Is state in atomic oxygen. 
Vertical line indicates the location of the oxygen nucleus. 

Figured shows the change in the electron density of the 
oxygen Is state in going from the oxygen atom to the wa- 
ter molecule. Again, the deformation of the Is states is 
quite clearly described in our multiresolution analysis de- 
spite being relatively small (0.16%, somewhat larger than 
in the oxygen molecule). In the water molecule, the elec- 
trons associated with the hydrogen atoms are stripped 
away and distributed in an approximately spherical shell 
of charge around the oxygen atom. Similarly to the wa- 
ter molecule, our calculation finds that this lessens the 
shielding of the protons and leads to the net attraction of 
the oxygen Is state toward each of the protons which is 
evident in the figure. The figure also shows that, in the 
water molecule, the displacement of the Is state leads 
to a net decrease of its amplitude on the nucleus. To 
quantify the impact of these core polarization effects on 
matrix elements describing transitions from the core to 
higher states, we have calculated the oxygen Is contri- 
bution to the dipole moment of the water molecule, find- 
ing 0.0003 clcctron-bohr. As a point of reference, the net 
dipole moment which wc compute for the entire molecule 
is 0.731 electron-bohr, in excellent agreement with the 
tabulated experimental value of 0.729 electron-bohr [ pO[ . 

Turning now to the energies associated with these 
states, we note that although common approximations 
to density-functional theory give relatively poor results 
for the absolute positions of deep core levels, our re- 
sults for the shifts in core-level energies show remarkable 
agreement with experiment. Accordingly, in comparing 
with experimental values, we maintain the relative place- 
ment of all of our eigenvalues (even between different 
molecules) by adding a single constant offset to our Is 
eigen-energies. Figure g| compares our multiresolution 
calculations of the energies of the oxygen Is levels with 
the experimentally observed X-ray photo-emission spec- 
trum of a mixture of water and molecular oxygen |p9| . 
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FIG. 8. Comparison of wavelet based all-electron calcu- 
lation of the location of oxygen Is states in the water and 
diatomic oxygen molecules (vertical bars) with experimental 
ESC A data (curve). Heights of bars for the O2 molecule are 
in proportion to the multiplicity of the spin-state resulting 
from the photo-emission. The height of the bar for H2O is 
arbitrary. ESCA data are after 

In molecular oxygen, the experiment shows significant 
splitting in the Is levels, with noticeably different cross- 
sections for the two levels. This results from the propaga- 
tion of paramagnetic bonding effects in the valence elec- 
trons down to the core levels. The stronger emission line 
results from the spin quartet state of the molecule (the re- 
sult of photo-emission from the minority spin channel), 
which has twice the expected cross-section of the spin 
doublet state (the result of photo-emission from the ma- 
jority spin channel). Our calculations, which treat the 
valence and core physics on a unified and equal foot- 
ing, predict precisely this valence-core effect. The figure 
displays our results for the oxygen Is Kohn-Sham eigen- 
values (with a single, constant shift for both molecules) 
as vertical bars with heights proportional to the cross- 
sections expected for each spin channel. For the oxygen 
molecule, we find precisely the correct splitting and or- 
dering for the two spin states. Finally, we note that our 
calculation of the relative positioning of the oxygen Is 
state in the water molecule is also in excellent agreement 
with the experiment. (As the relative strengths of the 
spectral lines between oxygen and water depend upon the 
relative concentrations in the experiment, the height of 
the bar in the figure indicating our calculation of the wa- 
ter molecule is arbitrary.) With this level of agreement, 
such calculations can be a great aid in the interpretation 
of spectra of greater complexity. 



VI. CONCLUSIONS 

Multiresolution all-electron calculations allow the uni- 
fied treatment of core and valence electrons where the 
local-density approximation is the sole approximation 
without systematic control. We have demonstrated that. 



with this approach, study of the impact of valence behav- 
iors on core electrons is straightforward and that mean- 
ingful core polarization effects and core-level shifts may 
be extracted to aid in the interpretation of experimental 
data such as ESCA. 

Further, we have shown that, with the use of newly 
developed techniques, such calculations may be carried 
out with an effort comparable to the corresponding pseu- 
dopotential calculations, thereby establishing for the first 
time multiresolution analysis as a viable approach to all- 
electron calculations. By accounting in detail for the 
computational demands of our calculations, we have es- 
tablished a standard against which the performance of 
future wavelet calculations should be compared. 
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TABLE HI. Restriction employed for all-electron calcula- 
tions of calcium: dimension (A^^, Ny,N^) and location relative 
to parent (A^,, Ay, Az) of each refinement level. 



level 
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N, 


Nz 
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TABLE IV. Restrictions employed for all-electron calcula- 
tion of oxygen molecule: dimension {Nx,Ny,Nz) and origin 
relative to its parent (A^,, Ay, A^) of each refinement level. 



level 


N, 


Ny 


A, 


A. 
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- 
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1 
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3 


1 


1 



TABLE I. Numerical representation of restricted multires- 
olution analysis in Figure hi: dimension {Nx,Ny) and location 
relative to parent (A^,, Ay) of each refinement level. 



a 
(a.u.) 



4V2Z 

Tz 



(a.u.; 



4-3.132576693428 Z 
-2.683558382240 Z 
-HQ.550981688812 Z 
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A, 
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TABLE V. Restrictions employed for all-electron calcula- 
tion of water molecule: dimension {Nx,Ny, N^) and origir 
relative to its parent (A^,, Ay, Az) of each refinement level. 



TABLE II. Widths a and norms Q of the three Gaus- 
sian charge distributions (5exp(— r-^/(2(T^))/(v'27ra)'' super- 
imposed to represent a nucleus of charge Z. All quantities 
reported in atomic units. 
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Quantity 


MRA 


V,s 


Nb 


200,000 


160,000 


Nf/b 


250 


1600 


Nf/t 


200- 10" 


ISO-IO" 


Na 


32,000 


8,500 


T 


8.2 hrs 


3.4 hrs 



TABLE VI. Analysis and comparison of computational 
demands of multiresolution analysis all-electron calculation 
(MRA) and plane-wave pseudopotential calculation (Vps) of 
oxygen molecule within LSDA for basis sets reaching a similar 
level of convergence. 
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